{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Translating Plate at 60°"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "using PotentialFlow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Set up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "At the beginning of every time-step, we first determine the bound vortex sheet strength required to satisfy the no-flow-through condition, then velocity of all vortex elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_ẋ! (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_ẋ!(ẋ, x, t)\n",
    "    plate, ambient_sys = x\n",
    "    motion, ambient_vel = ẋ\n",
    "    \n",
    "    Plates.enforce_no_flow_through!(plate, motion, ambient_sys, t)\n",
    "    \n",
    "    reset_velocity!(ẋ, x)\n",
    "    self_induce_velocity!(ẋ, x, t)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Once we have advected all the vortex elements, we release new segments of the vortex sheet from the edges of the plate.\n",
    "In this case, the strength of the segments are set to enforce the Kutta condition on both edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "shed_new_vorticity! (generic function with 3 methods)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function shed_new_vorticity!(sheet₊, sheet₋, plate, motion, t, lesp = 0.0, tesp = 0.0)\n",
    "    z₊ = (sheet₊.blobs[end].z + 2plate.zs[end])/3\n",
    "    z₋ = (sheet₋.blobs[end].z + 2plate.zs[1])/3\n",
    "    \n",
    "    segment₊ = Vortex.Blob.([sheet₊.blobs[end].z, z₊], [0.5, 0.5], δ)\n",
    "    segment₋ = Vortex.Blob.([sheet₋.blobs[end].z, z₋], [0.5, 0.5], δ)\n",
    "    Plates.enforce_no_flow_through!(plate, motion, (sheet₊, sheet₋), t)\n",
    "    \n",
    "    Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, segment₊, segment₋, t, lesp, tesp);\n",
    "    Vortex.Sheets.append_segment!(sheet₊, z₊, Γ₊)\n",
    "    Vortex.Sheets.append_segment!(sheet₋, z₋, Γ₋)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Discretize the plate and sheets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We are simulating a flat plate with zero thickness impulsively translating horizontally at a fixed velocity and angle of attack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "L = 2.0   # chord length\n",
    "N = 128   # number of plate control points (distributed along a extrema Chebyshev grid)\n",
    "\n",
    "ċ = 0.5L  # translation velocity\n",
    "α = π/3   # angle of attack\n",
    "\n",
    "Δt = 1e-2; # time step\n",
    "\n",
    "plate = Plate(N, L, zero(Complex128), α)\n",
    "motion = Plates.RigidBodyMotion(ċ, 0.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We place the initial segments of the vortex sheets perpendicularly off each edge of the plate.\n",
    "The sheets themselves are represented as discrete vortex blobs, all with the same blob radius `δ`.\n",
    "Both segments initially have unit circulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Vortex Sheet: L ≈ 0.01, Γ = 1.0, δ = 0.2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Δz₀ = im*Δt*exp(im*α)\n",
    "z₋, z₊ = plate.zs[[1,N]]\n",
    "\n",
    "δ = 0.2\n",
    "\n",
    "sheet₊ = Vortex.Sheet([0.5, 1.5] .* Δz₀ + z₊, [0.0, 1.0], δ)\n",
    "sheet₋ = Vortex.Sheet([0.5, 1.5] .* Δz₀ + z₋, [0.0, 1.0], δ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We then adjust the circulation of the vortex sheet segments to satisfy the Kutta condition.\n",
    "In this library, the vorticity flux from the edges of the plate is determined through the edge suction parameters.\n",
    "The Kutta condition simply corresponds to the suction parameters being zero at both edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Vortex Sheet: L ≈ 0.01, Γ = -0.489, δ = 0.2"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lesp = tesp = 0.0\n",
    "Plates.enforce_no_flow_through!(plate, motion, (), 0)\n",
    "\n",
    "Γ₊, Γ₋, _, _ = Plates.vorticity_flux(plate, sheet₊, sheet₋, 0, lesp, tesp);\n",
    "\n",
    "sheet₊ = Vortex.Sheet([0.5, 1.5] .* Δz₀ + z₊, [0.0, Γ₊], δ)\n",
    "sheet₋ = Vortex.Sheet([0.5, 1.5] .* Δz₀ + z₋, [0.0, Γ₋], δ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We organize the flow elements into the plate and ambient elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ambient_sys = (sheet₊, sheet₋)\n",
    "sys = (plate, ambient_sys)\n",
    "sys₊ = deepcopy(sys) # Used for storage during time-marching\n",
    "ẋs = (motion, allocate_velocity(ambient_sys));\n",
    "\n",
    "# Storage for pressure and forces\n",
    "pressures = Vector{Float64}[];\n",
    "Γs₋ = circulation(ambient_sys) .+ Plates.bound_circulation(plate);\n",
    "forces = Complex128[];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Time-Marching"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We use forward Euler to evolve the system and apply filtering on both the leading and trailing edge vortex sheets to surpress small-scale instabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "tf = 9\n",
    "T = 0:Δt:tf\n",
    "\n",
    "for t in T\n",
    "    plate, ambient_ω = sys\n",
    "    motion, ambient_u = ẋs\n",
    "    \n",
    "    TimeMarching.forward_euler!(sys₊, sys, t, Δt, compute_ẋ!, advect!, ẋs)\n",
    "    \n",
    "    # The force is computed based on the rate of change of impulse, which\n",
    "    # requires information about:\n",
    "    # - the plate's motion, \n",
    "    # - the strength, position, and velocity of the ambient vortex elements\n",
    "    # - the vorticity flux from the plate edges\n",
    "    push!(forces, Plates.force(plate, motion, ambient_ω, ambient_u,\n",
    "            (ambient_ω[1].blobs[end], ambient_ω[2].blobs[end]), Δt))\n",
    "    \n",
    "    # To calculate the pressure, we need the rate of change of the bound\n",
    "    # circulation distribution.  That also includes the bound vortex sheets\n",
    "    # generated by the newly shed vortex elements.  In order to account for\n",
    "    # that, we tell the `surface_pressure_inst` function where the new\n",
    "    # vortex elements are going to be released from so that it can compute\n",
    "    # their circulation.\n",
    "    push!(pressures, \n",
    "        let z₊ = (sys[2][1].blobs[end].z + 2sys[1].zs[end])/3\n",
    "            z₋ = (sys[2][1].blobs[end].z + 2sys[1].zs[1])/3\n",
    "    \n",
    "        Plates.surface_pressure_inst(plate, motion, ambient_ω,\n",
    "            (z₊, z₋), t, Δt, lesp, tesp)\n",
    "    end)\n",
    "    \n",
    "    sys₊, sys = sys, sys₊\n",
    "\n",
    "    # Redistribute the points along the free vortex sheet to be at least 3UΔt\n",
    "    # apart from each other.  At the same time, we filter out any length scales\n",
    "    # along the sheet shorter than 6UΔt\n",
    "    for s in 1:2\n",
    "        Sheets.arclength(sys[2][s]) ≥ 3Δt && Sheets.filter!(sys[2][s], 3Δt, 6Δt)\n",
    "    end\n",
    "    \n",
    "    shed_new_vorticity!(sys[2][1], sys[2][2], sys[1], ẋs[1], t + Δt, lesp, tesp)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Most flow elements in the library has [plot recipes](https://github.com/JuliaPlots/RecipesBase.jl) defined for them.\n",
    "So once [`Plots.jl`](https://github.com/JuliaPlots/Plots.jl) is installed along with the appropriate backend, we can directly call `plot` on the potential flow elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "using Plots\n",
    "clibrary(:colorbrewer)\n",
    "default(grid = false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "streamlines(-1:0.1:10, -3:0.1:2.5, sys, clim = (-3, 3), levels = collect(-3:0.3:3), color = cgrad([\"#efefef\", \"#efefef\"]))\n",
    "plot!(sys, color = :RdBu_r, clim = (-0.5, 0.5), legend = :none, cbar = :none, ratio = 1, colorbar_title = \"\\$\\\\gamma\\$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Forces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We computed the force on the plate in two different ways:\n",
    "\n",
    "- we computed the distribution of pressure differences across the plate along its control points.\n",
    "Since the plate control points are along an extrema Chebyshev grid, we can use the Clenshaw Curtis quadrature to integrate the pressure difference and find the normal force on the plate\n",
    "- we summed up the rate of change of impulse all the vortex elements in the system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let w = 0.5L*Plates.Chebyshev.clenshaw_curtis_weights(N),\n",
    "    labels = [\"Rate of Change of Impulse\" \"Integrated Pressure\"]\n",
    "    Fₚ = map(pressures) do p\n",
    "            -im*exp(im*plate.α)*(w[2:end-1] ⋅ p[2:end-1])\n",
    "        end\n",
    "    \n",
    "    plt = plot(layout = (2,1), size = (600, 400))\n",
    "    plot!(plt[1], T, 2*[imag.(forces) imag.(Fₚ)], \n",
    "          linewidth = 2, linestyle = [:solid :dash],\n",
    "          ylim = (0, 8), ylabel = \"\\$C_L\\$\", label = labels)\n",
    "    \n",
    "    plot!(plt[2], T, -2*[real.(forces) real.(Fₚ)], \n",
    "          linewidth = 2, linestyle = [:solid :dash],\n",
    "          ylim = (0, 10), ylabel = \"\\$C_D\\$\", xlabel = \"Convective Time\", legend = :none)\n",
    "    plt\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
